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In a recent paper [Phys. Rev. B 90, 115134 (2014)] we put forward a diagrammatic expansion for the self¬ 
energy which guarantees the positivity of the spectral function. In this work we extend the theory to the density 
response function. We write the generic diagram for the density-response spectrum as the sum of “partitions”. 
In a partition the original diagram is evaluated using time-ordered Green’s functions on the left-half of the dia¬ 
gram, antitime-ordered Green’s functions on the right-half of the diagram and lesser or greater Green's function 
gluing the two halves. As there exist more than one way to cut a diagram in two halves, to every diagram cor¬ 
responds more than one partition. We recognize that the most convenient diagrammatic objects for constructing 
a theory of positive spectra are the half-diagrams. Diagrammatic approximations obtained by summing the 
squares of half-diagrams do indeed correspond to a combination of partitions which, by construction, yield a 
positive spectrum. We develop the theory using bare Green’s functions and subsequently extend it to dressed 
Green’s functions. We further prove a connection between the positivity of the spectral function and the analytic 
properties of the polarizability. The general theory is illustrated with several examples and then applied to solve 
the long-standing problem of including vertex corrections without altering the positivity of the spectrum. In 
fact already the first-order vertex diagram, relevant to the study of gradient expansion, Friedel oscillations, etc., 
leads to spectra which are negative in certain frequency domain. We find that the simplest approximation to 
cure this deficiency is given by the sum of the zero-th order bubble diagram, the first-order vertex diagram and a 
partition of the second-order ladder diagram. We evaluate this approximation in the 3D homogeneous electron 
gas and show the positivity of the spectrum for all frequencies and densities. 

PACS numbers: 71.10.-w,31.15.A-,73.22.Dj 


I. INTRODUCTION 

Many-body perturbation theory (MBPT) has played an im¬ 
portant role in the understanding of the excitation properties 
of many-electron systems ranging from molecules to solids. 
An important class of excitations are the neutral excitations 
in which (in an approximate physical picture) electrons are 
excited from occupied to unoccupied states. These excita¬ 
tions can, for instance, be induced by external light fields 
and indeed the optical properties of materials, e.g., the in¬ 
dex of refraction, are completely determined by neutral ex¬ 
citations. For the understanding of the excitation spectrum 
many-body effects are of crucial importance as interactions 
lead to qualitatively new excited states of the system like 
plasmons and excitons in solids or the auto-ionizing states in 
molecules. In many-body theory the neutral excitation spec¬ 
trum is obtained from the density response function \ which 
can be calculated by diagrammatic methods. In practice one 
does not approximate \ directly but instead its irreducible part 
V, called the polarizability. The density response function 
and the polarizability are related through the integral equa¬ 
tion \ = V + Vvx where v represents the two-body inter¬ 
action between the electrons. The simplest approximation to 
V is the Random Phase Approximation (RPA) introduced by 
Bohm and Pines 1 to study plasmons in the electron gas. How¬ 
ever, many physical phenomena, such as excitons, are not de¬ 
scribed within the RPA. More complicated approximations in¬ 


volve typically a summation of an infinite class of diagrams, 
which is usually carried out with the Bethe-Salpeter equa¬ 
tion 2 . For instance, by summing ladder diagrams with a static 
screened interaction the excitonic properties of many solids 
are well described 2 . Nevertheless, there are several other cir¬ 
cumstances for which even the static ladder approximation is 
not enough. In fact, this approximation only allows for single 
excitations and therefore double and higher excitations are not 
incorporated 1 ' 4 . High excitations can significantly contribute 
to the spectrum of molecular systems. This is, for instance, 
the case in conjugated polymers, e.g., the polyenes, where the 
lowest lying singlet states have a double-excitation character 5 . 
Also for metallic systems there are several physical situations 
that require a theoretical treatment beyond the RPA 6 - 7 . The 
calculation of plasmon lifetimes 8 - 9 is one example. Another 
example is the calculation of the correlation induced double- 
plasmon excitations 6 ' 7 which has recently been studied both 
experimentally and theoretically for the simple metals. We 
also mention that beyond-RPA approximations have been in¬ 
vestigated actively in time-dependent density functional the¬ 
ory 10 (TDDFT), a widely used framework to study optical 
properties of molecules and solids. A key quantity within the 
TDDFT formalism is the so-called exchange-correlation ker¬ 
nel which depends entirely on corrections beyond the RPA. 
Several parametrizations of this kernel based on the homo¬ 
geneous electron gas exist 1 1-13 but their application to finite 
systems is problematic and the development of better approx- 
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imations is an active research area. 

From the viewpoint of diagrammatic MBPT the theoretical 
description of many-body interactions beyond RPA involves 
the inclusion of vertex corrections in the Feynman diagrams 
for the polarizability. Vertex corrections have been studied in 
several works 14-23 on the homogeneous electron gas. For this 
system it was found that first order vertex corrections give 
rise to negative spectral functions 17 ' 22 . Also for finite sys¬ 
tems it was found that by either restriction to a certain class 
of diagrams 3 - 24 or by truncation to a certain order in the iter¬ 
ation of the Hedin equations 25 the spectral function can be¬ 
come negative. This has two important consequences. First 
of all, it destroys the physical picture of the photo-absorption 
spectrum as a probability distribution and secondly it can lead 
to a density response function with poles off the real axis in 
the complex frequency plane thereby ruining its analytic and 
causal properties. Wrong analytic properties were, for in¬ 
stance, found in a study of atomic photo-absorption spectra 
with vertex corrections included to first order 24 . Concomi¬ 
tantly it was also shown that the absorption spectrum becomes 
negative around the energy of the inner-shell transitions. In 
this work we prove that a density response function with a 
positive spectrum has the correct analytic properties. 

The problem of negative spectra lies in the structure of the 
vertex correction and therefore the solution must be sought in 
the way diagrammatic theory is used. In a recent paper 26 we 
showed that a similar problem occurs in the spectrum of the 
Green’s function (which describes the photo-emission spec¬ 
trum rather than the photo-absorption spectrum). For that case 
we solved the problem by introducing the concept of “half¬ 
diagrams” which we then used to construct a diagrammatic 
expansion distinct from MBPT. In the present work we show 
that a similar procedure works for the polarizability too. We 
apply the theory to derive the lowest order vertex correction 
which preserves the positivity of the spectrum. We also evalu¬ 
ate this correction for the homogeneous electron gas by using 
a combination of analytical frequency integrations and numer¬ 
ical Monte-Carlo momentum integrations to evaluate the dia¬ 
grams. 

The paper is divided as follows. In Section IIA we intro¬ 
duce the spectral function for the density response function 
and the polarizability and give a diagrammatic proof of their 
positivity. The basic idea of the proof consists in cutting ev¬ 
ery diagram into two halves in all possible ways and in rec¬ 
ognizing the half-diagrams as the fundamental object of the 
diagrammatic expansion. Here we also derive the connection 
between the analytic properties and the positivity of the spec¬ 
tral function. In Section III we put forward a diagrammatic 
method to generate approximate polarizability with positive 
spectra. The method is first developed using bare Green’s 
functions and then extended to dressed Green’s functions. In 
Section IV we provide some illustrative examples of positive 
approximations and show how to turn a MBPT approxima¬ 
tion into a positive one by adding a minimal set of diagrams. 
We also address the positivity of approximations generated 
through the Bethe-Salpeter equation in Section V. In Sec¬ 
tion VI we evaluate the lowest-order vertex correction which 
yields a positive spectrum in the homogeneous electron gas. 


Our conclusions and outlooks are drawn in Section VII. 


II. THEORETICAL FRAMEWORK 


A. The density response function 


We study the properties of the reducible and irreducible re¬ 
sponse function within the Keldysh Green’s function theory. 
Although this theory is usually applied in the context of non¬ 
equilibrium physics we have found that it also provides a natu¬ 
ral and powerful framework for the calculation of equilibrium 
spectra. 

We consider a system of interacting fermions with Hamil¬ 
tonian 


H = J dx'i/^(x)/i(x)'0(x) 


+ - / dxdx' ft (x)$ (x')v(x, x , )^(x , )'0(x), (1) 


where the field operators ’ll;, ft annihilate and create a fermion 
at position-spin x = (r a), and v(x. x'j is the Coulomb in¬ 
teraction. The one-body part of the Hamiltonian is h(x) = 

— -ft- + qV (x), with V the scalar potential and q the fermion 
charge. Within the Keldysh formalism the correlators are de¬ 
fined on the time-loop contour C going from -oo to +oc 
(minus- branch C_) and back to — oo (plus- branch C+). Oper¬ 
ators on the minus-branch are ordered chronologically while 
operators on the plus-branch are ordered anti-chronologically. 
The Green’s function G(xi 0 i, x 2 z 2 ) (like any other two-time 
correlator) with times z\ and z 2 on the contour embodies 
four different functions G 0:3 , a. 3 = +/—, depending on 
the branch C_,C + to which z\ and z 2 belong. 26,27 For both 
times on the minus branch we have the time-ordered Green’s 
function G whereas for both times on the plus branch we 
have the anti-time-ordered Green’s function G ++ . The time- 
ordered and anti-time-ordered Green’s functions can be ex¬ 
pressed in terms of G~ + = G < and G + ~ = G > as follows 
(omitting the dependence on the position-spin variables) 


G ±± (G,f 2 ) = 0(ii - t2)G*(i 1 ,i 2 ) + 0(h - fi)G%,f 2 ). 

( 2 ) 

The four functions G° /7 form the building blocks of the fol¬ 
lowing diagrammatic analysis. 

The central object of this work is the density response func¬ 
tion x defined as the contour-ordered product of density devi¬ 
ation operators 

x(x 1 z 1 ,x 2 z 2 ) = -i(TclAn H (x 1 z 1 )An H (x 2 z 2 )]), 

where An(x) = n(x) — (n(x)), the subscript “ H ” implies 
operators in the Heisenberg picture and 7c is the contour¬ 
ordering operator. The average (...) is performed over the 
many-body state of the system. The greater x + ~ = X > anc l 
lesser x~ + = X < response functions read 

X > (xifi,x 2 f 2 ) = -i(Ah H (x 1 ti)An H (x 2 t 2 )) (3) 
X < (xHi,x 2 f 2 ) = -i(An H (x 2 t 2 )An H (x 1 t 1 )) (4) 
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and fulfill the symmetry relation 

*X^( x rU, x 2 ^) = [ix^(x-2t2,*iti)]* ■ (5) 


The quantity of interest for the excitation spectrum is the re¬ 
tarded response function 

X fl (xifi,x 2 f 2 ) = -i6(ti -f2)([nj?(xrfi),njT(x 2 t 2 )]) 

= 0(h - t 2 )(x > ~ X < )( x rfr,x 2 f 2 ). 

In equilibrium the Green’s functions G"^ as well as the 
response functions \“ /J depend on the time-difference t — i! 
and can therefore be conveniently Fourier transformed with 
respect to this time-difference. In the remainder of the paper 
we often suppress the dependence on spatial and spin indices 
and regard all quantities as matrices with one-particle labels. 
It is easy to verify that the Fourier transform of the retarded 
response function is given by 


f duj' B(uj') 

J 2tt lo — w' + irj 


with 77 a positive infinitesimal and spectral function 


( 6 ) 


%)=»[x > (w)-X < (w)]' (7) 


The matrix function B contains information on the energy of 
the neutral excitations and can be measured in optical absorp¬ 
tion experiments. From the relation in Eq. (5) we see that 
and ix < {oj) are self-adjoint and, therefore, B is self- 
adjoint too. The matrix function ix > {u>) vanishes for ui < 0 
and its average is proportional to the probability of absorb¬ 
ing light with frequency uj for ui > 0. In fact, *x > (w) is a 
positive semi-definite (PSD) matrix as it follows directly from 
the Lehmann representation. Similarly it can be shown that 
*X < (w) is PSD and vanishes for to > 0. Thus B(u>) is PSD 
for positive frequencies and negative semi-definite for nega¬ 
tive fequencies. Another property which follows directly from 
the definitions in Eqs. (3) and (4) is the relation ix < {u>) = 
[iX > (—w)]*. For Hamiltonians with time-reversal symme¬ 
try, such as in Eq. (1), this relation can also be written as 
iX < (oj) = [*x > (—w)] 1 ' which implies that B(ui) = —B(—oj). 
It is worth noticing for the present work that the PSD prop¬ 
erty is not guaranteed in diagrammatic approximations to the 
response function. How to construct PSD diagrammatic ap¬ 
proximations is the topic of the next section. 


B. Positivity of the exact response function 

The PSD property of the exact ix^(u) is manifest from 
the Lehmann representation of this quantity. It is instead 
less obvious to prove the PSD property from the diagram¬ 
matic expansion. Here we provide such a proof and bring 
to light a diagrammatic structure which forms the basis of a 
general scheme to construct PSD approximations (or to turn 
non PSD approximations into PSD ones by adding a minimal 
set of diagrams). We follow the same line of reasoning as in 
the recently published work on the PSD property of the self¬ 
energy . 26 The main difference is that the proof for the response 


function involves intermediate states consisting of particle- 
hole pairs rather than particle-hole pairs plus a particle or a 
hole. Since the derivation is otherwise similar we only outline 
the basic steps and refer to Ref. 26 for more details. Moreover 
since the derivation for x > is essentially the same as for x < 
we restrict the attention to \ < - 

The starting point is Eq. (4) for x < - Writing explicitly the 
time-evolution operator U in A fin we get 

*X<(1,2) 

= (T’o|W(fo,f 2 )An(x 2 )W(f 2 ,ii)An(xi)W(fi,f 0 )|^ , o)> 

where | H/o) is the non-degenerate ground state and the short¬ 
hand notation 1 = X| t\ and 2 = x 2 f 2 has been introduced. 
Under the adiabatic assumption the state (T'o) = U(to, r)|$o) 
is obtained by propagating the noninteracting ground-state 
|$o) from the distant future time r (eventually we take r —► 
00 ) to some finite time to with an adiabatically switched-on 
interaction. Therefore 

*X < ( 1 - 2 ) = ($o|W(r,f 2 )An(x 2 )W(t 2 ,T)^ 

i 

xW(r,fi)An(x 1 )W(f 1 ,r)|$ 0 ), (8) 

where we inserted the completeness relation in Fock space 
XT \<pi){<Pi\ = 1- The only states in Fock space which con¬ 
tribute in Eq. ( 8 ) have the form 

WpP) = c\ N -.c\ 1 c PN ..c Pl \^o), 

where dk, c\ annihilate and create a fermion in the fc-th eigen¬ 
state of the noninteracting problem. In Eq. ( 8 ) we can there¬ 
fore replace 

(9) 

i N—l pq_ 

where denotes integration or summation over the sets 

p = (pi,... ,pn) °f occupied states and q = (qi ,..., q n) 
of unoccupied states. The sum starts from N = 1 since 
the only state with N = 0 particle-hole pairs is |$o) and 
($o|Afi|$o) = 0. The prefactor in Eq. (9) originates from 
the inner product of the intermediate states 

= Sn > n ' 51 (-) P+Qs P(p)y S Q(q),q'’ 

PjQGttjv 

where P and Q run over the set 7 r..v of all possible permuta¬ 
tions of N indices and (—) p and (—)Q are the parities of the 
permutations P and Q respectively. Using Eq. (9) in Eq. ( 8 ) 
we can rewrite the lesser response function as 

00 

iX <(l,2)=^--^5W(2)5W*(l) ) (10) 

N=1 pq 

where the amplitudes S read 

= (^ ) |W(r,fi)Afi(xi)W(fi,T)|$o), (11) 
Spq\2) = ($o|W(t, f 2 )An(x 2 )W(f 2 , T )\ l Ppq ) )- (12) 
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The adiabatic assumption implies that turning the interaction factor: U(t, — t)|$o) = e IQ |$o). Hence we can rewrite the 
slowly on and off the state $0 changes at most by a phase amplitudes as (for more details see Ref. 26) 


SW*(xiit) 




($ 0 |T{e , f - dTH{T) cl 1 {T +).. ■ct jv (r+)c 9l (r). ■ ■ c qN (r) An(xifi)} |$ 0 ) 

($o|r{e-^-- d ^ ) }|$ 0 ) 

($ol T{e'^ d ^An(x 2 t 2 )cjjT) . ■ ■ 4(r)c PjV (r+).. ■ c Pl (r+)}|$ 0 ) 
($o|T{e l ^ d ^}|$ 0 ) 


(13a) 

(13b) 


with T and T the time-ordering and anti-time-ordering op¬ 
erators respectively. The time argument in the fermion cre¬ 
ation and annihilation operators specifies the position of the 
operators on the time axis, and r + denotes a time infinites¬ 
imally larger than r. Equations (13a) and (13b) show that 
S W* is an interacting time-ordered (N + 1)-Green’s func¬ 
tion whereas 5^1 is an interacting anti-time-ordered (AT+1)- 
Green’s function. Hence they can be expanded diagrammat- 
ically using Wick’s theorem. 28 The general structure of a S , 
S* diagram is illustrated in Fig. 1 and resembles half a \ dia¬ 
gram. The left-half corresponds to S^ N ’* with lines given by 
noninteracting time-ordered Green’s functions g whereas 
the right-half corresponds to with lines given by nonin¬ 
teracting anti-time-ordered Green’s functions g ++ (here and 
in the following we use the letter g to denote the noninteract¬ 
ing Green’s function). 

It is now easy to show that ix < (uj) is PSD. By Fourier 
transforming S with respect to t 2 and S* with respect to t-\ 
we find (omitting the dependence on the spatial and spin vari¬ 
ables) 



duj du' 
27r 27T 


x e -» t!+ »' tl ^ 5 w (w) 5 w, (w()i (14) 

pq 


In equilibrium \ < I s invariant under time translations, i.e., it 
depends on t\ — t 2 only. Imposing time translational invari- 



% 

N 



FIG. 1. (Color online) Typical diagram for S'*)].) (left) and S(2) 
(right) forming the lesser reducible response function. Labels p and 
q on the arrows denote quantum numbers of particles and holes re¬ 
spectively, see Eqs. (13). 


ance on the r.h.s. leads to 

OO 1 

E MM^ 5 ^ )(w)5 « 0 * (w ' )= - 77 ( w ) 5(a,-u/ ) (15) 

N= 1 ’ ’ pq 

with T some matrix function of the frequency to. Since for 
uj = oj' the l.h.s. is a sum of PSD matrices we conclude that T 
is PSD. Inserting Eq. (15) back into Eq. (14) we see that 
is the Fourier transform of i\ < which, therefore, is PSD too. 


C. Positivity of the exact polarizability 


In MBPT it is often more convenient to calculate the irre¬ 
ducible part V of x defined by the equation 

x{zi,z 2 ) = V{zi,z 2 ) 

+ J dzdz'V{z 1 ,z)v(z,z')x{z',z 2 ), (16) 

where all quantities are matrices in position-spin space and 
matrix product is implied. The two-body interaction in 
Keldysh space is given by v(z,z') = v5(z,z'). Diagram- 
matically V is obtained by removing from % all diagrams that 
can be separated into two pieces by cutting a single interac¬ 
tion line. Like the retarded response function in Eq. (6) the 
retarded polarizability has the spectral representation 



dw' H(w') 

27t uj — ui' + ig 


with spectral function 

B{u) =i[P>(uj) -V<(uj)}. 

The matrix functions V'z and (uj) are related by 


(17) 


(18) 


X $ H = [1 - V r (cj)v]- 1 V^(uj)[1 - vV a (uj)]-\ (19) 

where the advanced polarizability V a (uj) = \P r (uj)\* . As 
the exact is PSD this relation implies that is PSD too 
since [1 — V R v\ i = [1 — vV A ]. 

Alternatively we can use the diagrammatic expansion of the 
exact polarizability to show that V has the PSD property. In 
fact, the diagrammatic PSD proof for \ can easily be adapted 
for V. The removal of (interaction-line) reducible diagrams 
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from x is equivalent to the removal of 5-diagrams that can 
be separated into a piece containing the external ;\'-vertex (ei¬ 
ther 1 or 2) and a piece containing the pq vertices by cutting 
a single interaction line. We refer to these 5-diagrams as irre¬ 
ducible 5-diagrams or irreducible half-diagrams. It is easy to 
realize that the product between two irreducible half-diagrams 
yields an irreducible \ diagram. If we define 5 as the sum of 
irreducible half-diagrams then the lesser polarizability can be 
written as 

oo 

iV<(l,2)= £ MNiX g a (2)S £" (1) - (20) 

jV=l pq 

Fourier transfroming 5* with respect to t \ and 5 with respect 
to t- 2 , and carrying out the same analysis as in the case of y, 
see Eq. (14), we find that 

OO -j 

E (2D 

JV=1 " ‘ pq 

where is a PSD matrix function of the frequency ui. 

Since T is the Fourier transform of i'P < then iP < is PSD 
too. 

D. Partitions and cutting rule 

For the subsequent development of a PSD diagrammatic 
theory we need to introduce the concept of partitions of a V < - 
diagram. This concept naturally arises when we multiply two 
half diagrams, as we shall show below. Due to the anticom¬ 
muting nature of the fermionic operators we see from Eq. (13) 
that a permutation P of the p labels and a permutation Q of 

the q labels changes the sign of Spq 1 (and hence of 5^) by 
a factor (— ) P+ Q. Let then { Dpq } with j S Ijy be the set of 

topologically inequivalent diagrams for Sp^ that are not re¬ 
lated by a permutation of the p and q labels. By construction 
we have the expansion 

C = E E <22> 

j(zlN PiQ^N 

Inserting this expansion in Eq. (20) and using the fact that ttn 
is a group we get (for more details see Ref. 26) 

oo 

tp<(i,2) = jr £ £ r p+q 

JV=1 ji,j2 6/jv P,Q£kn 

xE^wO 1 )- (23) 

pq 

Next we observe that the in- and out-going Green’s functions 
with labels p and q are evaluated at the largest time r and, 
therefore, they are either a greater or a lesser Green’s function. 
At zero temperature the noninteracting lesser Green’s function 
g < satisfies the property 

E 3xpR T )Spy( r > ty) = ig*y(tx, ty) (24) 

V 


with a similar relation for the greater Green’s function. Hence, 
in the sum E pg D m\ 2 ) D pll P )Q 1 i yq) { 1 ) we can replace the 

product of two g£ with a single g^ connecting two internal 
vertices. The result is a polarizability diagram in which the 
lines of the left-half are time-ordered Green’s functions, the 
lines of the right-half are anti-time-ordered Green’s functions 
and the lines connecting the two halves are either lesser or 
greater Green’s functions. To represent this type of diagrams 
we label every internal vertex with a + or a — and introduce 
the graphical rule according to which a line connecting a ver¬ 
tex with label a = ± to a vertex with label /? = ± is a g a/3 . 
Let us name partition a V < diagram with decorated ± ver¬ 
tices. The full V < diagram is given by the sum of all parti¬ 
tions. The reverse operation of splitting a partition into two 
half-diagrams consists in cutting all the Green’s function lines 
between — and + vertices. In the following we refer to this 
reverse operation as the cutting rule. 

Before concluding this section we notice that in Keldsyh 
formalism a V < (ti,t 2 ) diagram is obtained from the corre¬ 
sponding Keldysh diagram V(z\, R by placing the contour 
time Z\ on the minus- branch C_, the contour time zq. on the 
plus- branch C+ and by integrating every internal contour time 
over C. Thus the V < diagram is the sum of diagrams with 
internal vertices decorated in all possible ways. Our deriva¬ 
tion shows that decorated diagrams which fall into multiple 
disjoint pieces by cutting the lines between — and + vertices 
do not contribute, i.e., they sum up to zero. In fact, these dia¬ 
grams cannot be written as the product of two half-diagrams. 


E. Connection between positivity and analytic properties 

It follows from the Lehmann representation that the re¬ 
tarded density response function x R ( UJ ) has poles at ± 0 , — ir] 
where f lj are the neutral excitation energies of the system. 
These poles lie just below the real axis in the complex fre¬ 
quency plane (in the continuous part of the spectrum these are 
smeared out to a branch cut). Therefore the function \ R {z) 
of the complex variable z is an analytic function in the up¬ 
per half plane Im z > 0. This analytic property is important 
in diagrammatic perturbation theory as the density response 
is an essential ingredient in the calculation of the screened 
interaction W as is, for instance, used in the GW approxima¬ 
tion. A violation of the analytic properties of \ R would lead 
to incorrect analytic properties of the Green’s function. It is 
therefore a relevant question to ask whether any approximate 
expression for the polarizability V R {uS) gives the correct ana¬ 
lytic properties of x R (u). We show here that this is the case 
whenever V R is PSD and the interaction matrix v is PSD as 
well (repulsive interaction). 

From the Dyson equation (16) for the response function we 
can write the retarded component in terms of the polarizability 
as 

X R (u)=V R (u)(l-vV R (cj)y 1 . (25) 

For our proof it turns out to be advantageous to look at the 
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function v 1 / 2 x R v 1 ^ 2 which can be written as 

v l/2 R v l/2 = V 1/2 V R ^)V 112 

X 1 - U 1 /27?fl( w ) w l/2 


where we used the PSD property of v to take the square 
root operation. The advantage of this expression is that 
v^V r {lu)v^ is PSD whereas vV R (u>) does not need to be 
PSD. Clearly since v 1 / 2 is frequency independent it does not 
change the analytic properties of x R (<*>)■ 

We want to know whether X R ( Z ) has a pole in upper half 
of the complex plane, i.e., for z = x + iy with y > 0. This 
is either the case when V R (z) has a pole at z or when 1 — 
v l / 2r p R (z)v i/ 2 has a zero. From Eq. (17) we have 



duj' B(u') 

2n x — to' + iy 


(27) 


and therefore 


Im V R (z) = -y 


R eP R {z) 



f dco' B(u') 

J 2tt (x — w') 2 + y 

dco' (x — uj')B(oj') 

2tt (x — to') 2 + y 2 ’ 


2 > 


(28) 

(29) 


where we used the short-hand notation lmV R = (V R — 
V a )/ 2i and ReV R = ( V R + V A )/2. Since y > 0 and B 
is an integrable function both integrals are finite and V R (z) 
does not have a pole. It remains to consider the possibility 
that the operator 1 — v 1 ^ 2 'P R (z)v 1 ^ 2 has a zero eigenvalue. If 
this was the case then there would exist an eigenvector |A) for 
which (l — v 1 ^ 2 V R {z)v 1 ^ 2 ) | A) =0, and hence 

(A| 1 - v 1/2 ReV R {z)v 1/2 |A) = 0, (30a) 

(A| v 1/2 lmV R {z)v 1/2 |A) = 0. (30b) 

Using Eq. (28) we rewrite Eq. (30b) as 
{X\v 1 / 2 hnV R (x + iy)v 1 / 2 \X) 

= -yj~ ^<A|t(31) 


in contradiction with the assumption. We conclude that x R { z ) 
cannot have poles in the upper half plane when x = Re 2 ^ 0. 
This leaves us with the case x = 0. For x = 0 the function 
£(ui) = 0 and Eq. (30b) is automatically satisfied. Instead 
Eq. (30a) reads 


(A 11 — v 1/2 Re'P' R (z)u 1 / 2 |A) 

f 00 ^/ (A|u 1 / 2 ff(u;> 1 / 2 |A)cy 
+ Jo 2 tt u}' 2 + y 2 


(33) 


since v Y l 2 Bv x l 2 is PSD for positive frequencies. In this case 
too Eq. (30a) cannot be satisfied and therefore x R cannot have 
poles in the upper half plane when V R and v are PSD. 

We mention that the situation is different for negative 
semidefinite interactions. In this case a similar formula can be 
defined but with a minus sign in front of the integral. In fact, 
for a strong enough attraction the right hand side of Eq. (33) 
can be zero for a well-chosen value of y and consequently 
X R can have poles at z = iy in the upper half plane. The 
violation of the analytic property occurs for instance in the at¬ 
tractive Hubbard dimer 29 . As a final remark we note that in a 
similar fashion one can prove that when the spectral function 
of the self-energy is PSD then the Green’s function is analytic 
in the upper half of the complex frequency plane. 

Implications to the f-Sum Rule. The connection between 
positivity and the analytic structure has important conse¬ 
quences for the f-sum rule, which relates the first momentum 
of the retarded density response function to the equilibrium 
density no( x ) 


Lx fi (x,x'u) dut = -?7rV' • V (n 0 (x')S(r - r')) 

(34) 

The derivation of the /-sum rule assumes that the integral over 
the retarded response function can be closed on the upper half¬ 
plane 27 , i.e., x R is analytic in this region. In accordance with 
the results of this Section the analytic assumption is verified 
for those approximations which fulfill the PSD property. 


III. PSD DIAGRAMMATIC EXPANSION 


where we took into account that B(ui) = —B{—u>) to write 
the integral between 0 and oo, and defined the function £{ui) 
according to 


(x — oj) 2 +y 2 {x + w) 2 + y 2 \ 

The function £(u) is odd and positive (negative) for posi¬ 
tive frequencies and for positive (negative) x values. There¬ 
fore £(cu) has a definite sign on the positive frequency 

axis when x is nonzero. Let us first consider this case, 

i.e., x ^ 0. Since we assumed that y > 0 and since 

(A|u 1 ^ 2 H(w , )u 1 / 2 |A) > 0 for positive frequencies the only 
way to have (A| v 1 ^ 2 lmV R (z)v 1 ^ 2 |A) = 0 is to demand that 
(A|u 1 / 2 H(aj , )'t; 1 / 2 |A) = 0 for all u/. This however, would 
imply from Eq. (29) that also (A|u 1 / 2 Re7 : ’' R (z)u 1 / 2 |A) = 0 
which in turn would imply that (A| 1 — v 1 / 2 V R {z)v 1 / 2 1 A) = 1 


A. Formulation with noninteracting Green’s functions 

In Section IIA we have given a diagrammatic proof of the 
PSD property of the exact polarizability. We further showed 
that this PSD property is essential to guarantee the correct an¬ 
alytic structure of the density response function. In practice, 
however, the polarizability is calculated in a MBPT fashion 
by considering a subclass of Feynman diagrams and hence 
the PSD property is not, in general, a built-in property of the 
approximation. The proof given in Section IIC paves the way 
for a diagrammatic theory of PSD polarizabilities which is al¬ 
ternative to the more standard MBPT. We have seen that the 
PSD property follows from the formation of perfect squares 
and that these squares are the sum of partitions. In contrast 
MBPT gives a sum of diagrams where each diagram is a sum 
of partitions but, in general, do not form perfect squares. The 
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most general MBPT approximation to the polarizability when 
written in terms of partitions (or equivalently in terms of prod¬ 
ucts of half-diagrams) reads 

OO 

»P<(1,2 )=X; £ £ R P+Q 

N—i pc n UiJ2) 

Q&'SU* 

xE^ a) ( 2 )ft (!) ( 1). (35) 

where Xjv is a subset of the product set / y x /.y and for any 
given couple (ji, j 2 ) the sums over P and Q run over a subset 
TT^'p 2> and Tr^'J 2 " 1 of the permutation group 7 Tjv. The mini¬ 
mal number of additional partitions to add in order to turn the 
MBPT approximation into a PSD approximation is found by 
imposing on T >< the mathematical structure in Eq. (23). This 
is achieved as follows, see also Fig. 2. Let {/^} be a set of 
disjoint subsets of / y with the property that the union of the 
product sets 

\Jl%xI%Dl N (36) 

a 

and contains the least number of elements of In x / y. Since 
the subsets Jjy are disjoint the sets 1^ = In IT (Ift x In) are 
disjoint too and due to Eq. (36) we have that |J Q Xjy = ^-n- 
For any given a we then consider the smallest subgroups of 
the permutation group ttn with the property that 

*N,p T U ^n, p 32 \ (37) 

(jih)£Z% 

R, => U ( 38 > 

By construction the polarizability 

OO 

^s»(u) = EE E E <-) p+0 

N=1 a 3ij2ei% p £*N, P 
Q^N,q 

ED«'>(2)C^ ( „(l) (39) 

pq 

contains all partitions of Eq. (35) plus the minimal number 
of additional partitions to form perfect squares. Consequently 
the polarizability in Eq. (39) is a diagrammatic PSD approxi¬ 
mation. More precisely we can say that a PSD diagrammatic 
approximation to V is not the sum of MBPT diagrams, rather 
it is the sum of partial or decorated diagrams (the partitions) 
with internal vertices either on the minus or the plus branch 
of the Keldysh contour. Examples of the PSD procedure are 
presented in Section IV. 

As a final remark we mention that another cutting proce¬ 
dure based on time-ordered half diagrams has been used by 
Sangalli et air for the Bethe-Salpeter kernel of the so-called 
second RPA approximation. In that work the Lehmann prod¬ 
uct structure was important to satisfy an identity for the de¬ 
terminant of the Bethe-Salpeter kernel which guarantees that 



FIG. 2. (Color online) Decomposition |J a Xj(r = In of the In 
subset (denoted as three disjoint gray areas) into a union of (two in 
this example) product sets. 

no spurious poles occur in the density response function ob¬ 
tained by solving the Bethe-Salpeter equation. However, the 
procedure proposed in Ref. 3 cannot be used to obtain a direct 
expression for the spectral function, thus making it difficult 
to address the PSD property. The difficulty has its origin in 
the fact that the standard time-ordered formalism is not the 
natural formalism to express the spectral function. As it has 
been shown in this work the Keldysh contour technique facil¬ 
itates enormously the calculation of the spectral function of 
any diagram, the result being a product of time-ordered and 
anti-time-ordered half diagrams joined by lesser and greater 
Green’s function lines. 


B. Formulation with dressed Green’s functions 

In the previous Section we discussed how to generate PSD 
diagrammatic approximations for the polarizability. A dia¬ 
grammatic approximation is the sum of partitions and for the 
PSD property to be satisfied the product of half-diagrams re¬ 
sulting from the cut partitions has to form the sum of per¬ 
fect squares. The possibility of cutting a partition relies on 
Eq. (24), according to which the product of two g> yields a 
single g>. This property, however, is valid only for noninter¬ 
acting Green’s functions. It would be extremely useful to for¬ 
mulate a cutting rule for partitions written in terms of dressed 
Green’s functions. The main advantage of a dressed PSD for¬ 
mulation is the absence of polarizability-diagrams with self¬ 
energy insertions, thus enabling us to work exclusively with 
skeleton diagrams. 

Here we consider partitions with dressed Green’s function 
lines and show how to write these partitions as the product of 
two half-diagrams. We therefore need to replace Eq. (24) with 
some other equation where G^- is expressed as the “product” 
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of two functions. The main idea has been discussed in detail 
in our recent work on the self-energy. 26 In frequency space the 
greater and lesser Green’s function read 

G^(xifi,x 2 f 2 ) = =R J ^A^(xi,x 2 ;w)e - *“ (tl ~ t2) , 

(40) 

where A < (oj) = f(ut)A(u}) is the removal part of the spectral 
function A(uj) whereas A > (uj ) = (1 — f(oj))A(i>j) is the ad¬ 
dition part of A(uj), and f(u) is the zero temperature Fermi 
function. If the self-energy is PSD then both A > and A < are 
PSD. We expand the matrix A ^ (w) in terms of its eigenvalues 
an (uj) >0 and eigenvectors u n (uj,x) 

A^ (xi, x 2 ; oj) = ^a|(w)M n (w,Xi)w*(w,x 2 ) (41) 

n 

and define the square root matrix as 



n 


We then make the rule that when cutting a partition the in¬ 
ternal lines are G~ for the left half, G ++ for the right half, 
and 


V / G^(xiG,x 2 f 2 ) = /^(xi,x 2 ;.)e lul{tl * 2) 

(43) 

for the dangling lines of the two halves. For the reverse op¬ 
eration of gluing the two halves we make the rule that the 
“product” of two dangling lines is defined according to 

J dydt v / G< Xiy (G,i)v / G< y x 2 (f,f 2 ) =iG < ( 1 ,2), (44a) 
JdydtVG> Xiy (t 1 ,t)VG> yX2 (t,t 2 ) = -iG > { 1 ,2). (44b) 

These equations replace Eq. (24) and the analogous for g > 
in the dressed case. Except that for an additional time- 
integration Eqs. (44) have the same matrix-product structure 
as in the undressed case and it is straightforward to verify that 
the gluing of two half-diagrams gives back the original par¬ 
tition. This implies that if a certain sum of partitions with 
undressed Green’s functions g is PSD, and hence it can be 
written as in Eq. (39), then the same sum of partitions with 
dressed Green’s functions G is PSD too since the correspond¬ 
ing polarizability can be written as 


iV, 


-m-EE E E <-) 


P+Q 


iV=l a 


pz* n, p 

Q^N.g 


x / dl„ ... dN„ / dl n ... dN„ 


X D lp.’..N p ,l q ...N q ( 2 ) D ptl p ...N p )Q(l q ...N q )( 1 )’ ( 45 ) 


where we introduced the short-hand notation n p = ( p n , t n ' 1 ) 

(n = 1 ,N) and m q = ( q m , t[n ) (m = 1,..., N) as well 


o 



FIG. 3. (Color online) Partition and decomposition in half-diagrams 
of the RPA bubble diagram. 


as 

dt$. (46) 

Pn " " 9m ” 

In Eq. (45) the D’s represent the half diagrams with dangling 
lines Vg $ and external vertices in l p ... N p , l q ... N q . 


dn p = ^2 / dt n ) I / dm q = ^2 


IV. EXAMPLES 

In this section we consider some commonly used diagram¬ 
matic approximations to the polarizability and address the 
PSD property for each of them. 

a. Zeroth order In Fig. 3 we show the RPA bubble. This 
diagram can be partitioned in only one way and the decom¬ 
position in terms of half-diagrams is shown on the right of 
the equality sign. According to the cutting rules for dressed 
diagrams we have 

D\ plq { 1) = V / G< XlPl (G, 4 p) )\/G>g lXl (4 9) , G), 
Di p i q ( 2) = VG< P1X2 (4 p) , t 2 ) VG> X2qi (t 2 ,t[ q) ), 

and it is straightforward to verify that the RPA bubble can be 
written as 

7>< a (1,2 )=~ij d l p dl q D lplq (2)D* lplq (l). (47) 

The RPA V has clearly the structure of Eq. (45) and therefore 
the spectrum for the response function is positive. 

b. Simplest vertex In Fig. 4(a) we consider the lowest 
order (in bare Coulomb interaction) vertex diagram for V. 
The cutting rules yield only two partitions (on the right of the 
equality sign) since the bare interaction is local in time and 
hence = 0. From the decomposition in half-diagrams 
we see that the vertex diagram does not have the structure 
of Eq. (45). Therefore the PSD property is not guaranteed. 
The minimal set of partitions to add follows directly from the 
rules of Section III A. The result is shown in Fig. 4(b). By 
multipling the half-diagrams in the curly brackets we get two 
additional diagrams: the RPA bubble and a partition of the 
second-order ladder diagram. Here and in the following we 
use the convention that an internal vertex without —/+ labels 
implies a summation over —/+. From this example we also 
infer that the sum of only the RPA bubble and the vertex dia¬ 
gram is not, in general, PSD. Noteworthy this sum constitute 
the so called exact exchange (EXX) approximation to the ker¬ 
nel of Time-Dependent Density Functional Theory (TDDFT). 
The EXX kernel has been calculated in Ref. 24 for the case 
of closed shell atoms and it was found that it has poles in 
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the upper half of the complex frequency plane. This incor¬ 
rect analytic behavior is a direct consequence of the relation 
between the PSD and the analytic properties, as discussed in 
Section IIE. The correct analytic properties of the EXX ker¬ 
nel can be restored by adding the partition of the second-order 
ladder diagram shown in Fig. 4(b). 

c. RPA screening In bulk systems the electron screen¬ 
ing plays a crucial role and, therefore, one usually works with 
interaction-skeletonic diagrams in which the bare interaction 
v is replaced by some screened interaction W. Below we ap¬ 
ply the theory developed in Section III A to a few diagram¬ 
matic approximations in which W is the RPA screened in¬ 
teraction, i.e., W = v + vV Q W where Vq = — iGG . Tak¬ 
ing into account that partitions with isolated +/— islands do 
not contribute to the polarizability we immediately conclude 
that W and W ++ should not be cut along the internal 
Green’s function lines. For the lesser/greater H’-lines we note 
that W< = W VqW ++ and W> = W++V^W—, see 
Fig. 5. Therefore the cut of a amounts to a cut of the two 
Green’s function lines in Vq". 

The simplest diagrammatic approximation in terms of W 
is given by the four partitions of Fig. 4(b) in which v —» W, 
see Fig. 6. We observe that the sum of them does not give the 
full vertex diagram since the partitions with W h = W < and 
W+~ = W> are missing. These partitions would vanish if, 
instead of the RPA W, we used an externally given static W 
like, e.g., a Yukawa-type interaction. We will come back to 
this example in Section VI where we give explicit expressions 
of the diagrams of Fig. 6 for the electron gas and show the 
crucial role played by the second-order ladder diagram for the 
PSD property. 

cl. First order The decomposition into half-diagrams of 
the full first order (in screened Coulomb interaction) vertex 
diagram is shown in Fig. 7(a). Since the times of the internal 



FIG. 4. (Color online) (a) Partition and decomposition in half¬ 
diagrams of the vertex diagram. The wiggly line denotes the bare 
interaction v. (b) The minimal set of additional partitions to restore 
the PSD property. 



FIG. 5. Decomposition of the screened interaction (thick wiggly line) 
W < in half-diagrams. The decomposition of W > is analogous and 
given by the bottom diagram with — +. 

vertices can be different (nonlocal interaction) we have four 
different partitions. In two of them we only cut the G-lines 
while in the other two we also cut the IT-line (the cut of the 
W- line leads to half-diagrams with four dangling lines). Thus, 
we have two half-diagrams with one particle-hole (D^ and 
D^) and two half-diagrams with two particle holes ( D ^ and 
D ( & )). Naming each half-diagram as shown in the figure 7(b) 
we can write (omitting the integrals over the vertices to be 
glued as well as dependence on the external vertices) 

iV< = £>(°) £>W’ + Z)( b ) £>(“)* + D^)f)^ b y + D w D {ar , 

(48) 

which does not have the structure of Eq. (45) and hence it 
is not PSD. Applying the rules of Section III A we find the 
minimal set of partitions to add in order to restore the PSD 
property 

^PSD = E £>(■?)*+ D {i) D U) \ (49) 

ij=a,b ij=a,b 

where in each of the sums the indices i and j independently 
take values a, b. The resulting diagrammatic approximation 
to V is illustrated in Fig. 7(c). The important message of 
this example is that the additional diagrams are not necessar¬ 
ily skeletonic in G (occurrence of self-energy insertions, viz. 
the last two diagrams in the figure). Thus attention has to be 
paid when restoring the PSD property using a dressed G. In 
order to avoid double countings one should not dress the G 
with the same self-energy appearing in the diagrams of the 
PSD polarizability. 

e. GW exchange-correlation kernel We conclude this 
section with another important example. In Ref. 30 it was 
shown how to generate conserving approximations to the 
TDDFT kernel using a variational principle a la Luttinger- 

o + vj> + + xtE 

+ + 

FIG. 6. (Color online) The simplest PSD approximation to the polar¬ 
izability in terms of the screened interaction W (thick wiggly line). 
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FIG. 7. (Color online) (a) Decomposition of the vertex diagram with 
screened interaction line into half-diagrams, (b) Constituent half¬ 
diagrams. (c) The resulting PSD polarizability. 


FIG. 8. (Color online) (a) Diagrammatic approximation to the polar¬ 
izability. (b) Constituent half-diagrams in addition to those presented 
in Fig. 7(b). 


half-diagrams with one particle-hole and three half-diagrams 
with two particle-holes. Some of them have already been in¬ 
troduced in Fig. 7(b), the new ones are shown in Fig. 8(b) 
and the expression of V < in terms of them is (again omitting 
integrals and the dependence on the external vertices) 

iV< = [b { p ] Q + D W + D% + D { p ] + D * 

+ D W ( b { p Q + DM + D ( p ] Q + D ( p ] + * 

+ D (°) (. D& + £>pq + D { p ] + D ( q ] + D { p ] + D * 

+ D (a) (l) (fe) + D {c) + D {d) y 

+ + D (c) + L> (a) *, (50) 

where the half-diagrams with subindex P (and/or Q) are cal¬ 
culated at permuted values 2 p , l p (and/or 2 q , l q ). For in¬ 
stance, D^DpQ is a short form of the integral 


Ward. 31 The underlying variational functional of Luttinger- 
Ward is a functional of the bare interaction v and the Green’s 
function G. To lowest order in v one can show that the vari¬ 
ational principle leads to the EXX approximation. It is pos¬ 
sible to extend the Luttinger-Ward idea to functionals of the 
screend interaction IT' and the Green’s function G. 31 In this 
case the lowest order approximation is the “time-dependent 
G IT” (TD GW) approximation; and the TDDFT kernel (more 
precisely its convolution with two 'Pq’s) is given by the sum 
of the diagrams in Fig. 8(a) (see also Sec. III.B of Ref. 30). 
These are the same diagrams evaluated by Stememann et al. 1 
and by Huotari et al. for the electron gas in order to explain 
the double-plasmon shoulder in the absorption spectrum of 
sodium. By partitioning each diagram of the approximate po¬ 
larizability we find that V < can be written in terms of four 




dl q d2 q d[^2 p 1,2, ( 2 ) -02 d 1 d 2 o 1„( 1 )- 


4p±p*q± q 


The polarizability of Eq. (50) is not, in general, PSD since 
it does not have the form of Eq. (45). The minimal addition 
to restore the PSD property follows from the rules of Sec¬ 
tion III A which give 


iPpsD= E dWdW+ £ £ D^D%. 

ij=a,b,c,d ij=a,b,c P,Q£tv 2 

(51) 

The first sum leads to 4 x 4 = 16 partitions whereas the second 
sum leads to 3 x (2 2 x 3) = 36 partitions. In Fig. 9 some 
representative ones are shown. 
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FIG. 9. (Color online) A few additional partitions of Eq. (51). 


V. ON THE POSITIVITY OF THE BETHE-SALPETER 
KERNEL 

So far we have studied only approximations to the irre¬ 
ducible response function consisting of a finite number of di¬ 
agrams. However, typically approximations beyond the RPA 


involve an infinite series of diagrams conveniently resummed 
through the Bethe-Salpeter equation (BSE) 2 - 33 . A natural 
question to ask is whether the corresponding spectral function 
is PSD. To answer we have to find the diagrammatic structure 
encoded in the kernel of the BSE. The BSE is a Dyson-like 
equation for the four-point reducible polarizability £(12; 34) 
and it is obtained as the response to a non-local scalar poten¬ 
tial u( 4, 3) 2 

£(i2;34 >-^l= £ » (i2!34) 

+J d(5678)G(l, 5)G(7,3)/C(56; 78)£(82; 64) (52) 

where £ 0 (12;34) = G(l, 4)G(2,3) and the four- 

point reducible kernel 1C is given by /C( 12; 34) = 

—*5£(1,3)/5G(4,2). The variation of the Green’s function 
is related to the two-particle Green’s function and therefore C 
is related to the two-particle excitation spectrum. By taking 
the limit 3 —> 1 + and 4 —>• 2 + we obtain an equation for the 
response function since x(l, 2) = iC( 12; l + 2 + ). 

From standard approximations to the self-energy, e.g., 
Hartree-Fock (HF), second order Born (2B) or the GW ap¬ 
proximation, we can derive a diagrammatic expression for 
the kernel 1C, see Fig. 10. By defining a two-particle irre¬ 
ducible and one interaction line irreducible kernel 1C as (see 
also Fig. 11(a)) 

£(12; 34) = /C(12; 34) - iS( 1, 3)5(2,4)v(l, 2), (53) 



JC 


HF 


: Vwv4 + 


we can write the polarizability as an infinite series of response 
diagrams 

V{1,2 ) = p 0 (l, 2) + ij d(3456) G(5,1)G(1,3)£(34; 56) 


(b) o 

^2B ~ | + 




/C 2B = + I + 


| o | * o * M 


+ 




(c) 9 

^GW = \ + + 



^GW = + 


+ 




xG(4,2)G(2,6) + ... (54) 


where ^(l, 2) = — *G(1, 2)G(2,1). The diagrammatic ex¬ 
pression for this equation is shown in Fig. 11(b). By using 
the kernels in Fig. 10 we obtain the approximations shown in 
Fig. 12. The HF approximation for the BSE kernel yields the 
diagrams shown in the Fig. 13, which we can easily see to 



FIG. 10. (Color online) Self-energy and the corresponding BSE 
kernel for (a) Hartree-Fock (b) Second order Born and (c) GW ap¬ 
proximations. 


FIG. 11. (Color online) (a) Two-particle irreducible and one interac¬ 
tion line irreducible kernel 1C. (b) The diagrammatic expression for 
the polarizability in terms of BSE kernel 1C. 
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< b > p -=o + <X> + <D> 


FIG. 12. (Color online) Approximations for the polarizability by 
using BSE kernel with various self-energy approximations, (a) 2B 
polarizability, (b) GW polarizability. 


be PSD. From this example we also deduce that if a kernel 
cannot be partitioned then the corresponding polarizability is 
PSD. A commonly used approximation to study the exitonic 
properties of solids is the static GW approximation— 3,29 with 
kernel 

4w(12; 34) = iS( 1,2)<5(3, 4)W(l, 3), (55) 

where the functional derivative of W with respect to G is 
neglected. This approximation leads to the polarizability of 
Fig. 13 in which the bare interaction lines are replaced by 
statically screened ones. Therefore, the static GW approxi¬ 
mation yields PSD spectra. For the full GW approximation 
some of the half-diagrams have been already worked out in 
Fig. 8 and the resulting diagrams after the gluing procedure 
have been shown in Fig. 9. From this figure we see that the 
PSD procedure leads to new types of diagrams which are not 
obtained via the iteration of the BSE. Thus, we conclude that 
the BSE polarizability with GW kernel does not necessarily 
have PSD spectra. The same is true for the 2B approxima¬ 
tion as can be seen from the half-diagrams generated by the 
last diagram on the first line of Fig. 12. By applying the PSD 
procedure we will generate diagrams which are not obtained 
via iteration of the BSE. Therefore, even the 2B approxima¬ 
tion is not a PSD approximation for the BSE. These simple 
examples show that the kernels generated by conserving «1>- 
derivable self-energies 31,34, 0 do not need to be PSD. 


VI. NUMERICAL RESULTS 

As an illustration of our method we compute the spectral 
functions B{k,ui) of the polarizability V(k,u) for the three- 
dimensional homogeneous electron gas. For convenience we 
introduce here the spectral functions for the positive B > (k, oj) 
and negative B < (k, oj) frequencies and scale them by the fac¬ 
tor of 8nar s in order to make the zeroth order (the Lindhard 


p h F = o + cd + <y> 



FIG. 13. (Color online) The polarizability calculated from the HF 
BSE kernel is PSD. 


polarization function) density independent 


k, 
k, 

(56) 

where we expressed k in units of the Fermi momentum ky 
1 / (ar s ) and u> in units of the Fermi energy ty = k\ /2. Here 
a = [4/(97r)] 3 / 3 and r s is the standard measure of the sys¬ 
tem’s density-the Wigner-Seitz radius-expressed in units of 
the Bohr radius. 

In 1958 J. Hubbard diagrammatically studied the correla¬ 
tion energy of a free-electron gas 36 and introduced what is 
now known as the local field factor (f(k,u>), same notation 
as in the original manuscript is used). A very interesting in¬ 
troduction to the historical development of this concept and 


B {0) (k,oj ) = 


-4 


M <2- 

k — 


\2-k\ < M <2 + 



k / k F 


FIG. 14. (Color online) Distribution of positive (green) and negative 
(pink) values of the first order spectral function of the vertex diagram 
in Fig. 4(b) in the k — ui plane. The lines represent isospectral curves 
with values ±0.3 (dashed) ±0.1 (dotted) and 0 (solid). 
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FIG. 15. (Color online) Sum of these zero to second order polariz¬ 
ability diagrams yields positive spectral function. Lines with arrows 
denote the electron propagator Go(k,co), whereas wavy-lines stand 
for bare or screened Coulomb interaction. Vertices are labeled with 
+ (—) if they belong to the positive (negative) time ordering part of 
the Keldysh contour. 


its importance for the density functional theory can be found 
in Ref. 37. In essence, it provides a simple way to go beyond 
RPA both in the treatment of the density response function and 
the total energy of a many-body system and relates between 
the exact and zeroth order polarizabilities v(k)f(k,co) = 
[V(k, cu)]- 1 - [P<-°)(k, cu)]- 1 . Thus, every advancement in 
the calculation of the proper density response function leads 
to our improved knowledge of the local field factor and, corre¬ 
spondingly, of the density functionals. After Hubbard’s origi¬ 
nal static approximation f(k,co) ~ k 2 /(k 2 + k 2 F ) there were 
numerous works to compute the local field factor using the 
above relation, notably the exact long and short wavelength 
limits (see Ref. 37 and references therein). Diagrammatically, 
the simplest case of the first order diagrams (in terms of the 
bare Coulomb interaction) was computed in a concise form by 
Engel and Vosko 23 for the static case. Importantly, they con¬ 
sidered two types of diagrams, the proper first order response 
and the first order self-energy insertions and demonstrated a 
rather large cancellations between these two contributions. In 
a full generality the frequency dependent first-order results 
were obtained by Holas, Aravind, and Singwi 17 . However, 
the analytic form is rather complicated and can only be ex¬ 
pressed in terms of a one-dimensional integral. We will use 
these results for a comparison and, therefore, numerical de¬ 
tails concerning the evaluation of this integral are presented 
in Appendix A. In fact, already the first order vertex diagram 


of Fig. 4(b) demonstrates the problem of standard MBPT. In 
Fig. 14 we depict its momentum and energy resolved spec¬ 
tral function B^ computed according to the expression of 
Ref. 17. The pink shaded area denotes a part of the particle- 
hole continuum where the 1 st order spectral function is nega¬ 
tive. 

In order to solve this problem we use our method and 
consider the irreducible polarizability diagrams shown in 
Fig. 4(b) and Fig. 6: a particular second order diagram must 
be added in order to compensate for the negative sign of B^\ 
The two sets of diagrams are topologically identical, however, 
the first one is given in terms of bare Coulomb, while the sec¬ 
ond contains screened interacting lines. To be general we will 
start with the second more complicated case and derive the 
first case by making the limit w(k) —> fl(fc) —► oo in the 
plasmon pole approximation for the screened Coulomb inter¬ 
action: 


W 0 (fc,w) 


v(k) w(k ) w(k) 

2 co — fl(fc) + irj co + f l(k) — irj 


(57) 

In this expression w(k) = f(fc)fl 2 (0)/f2(/c), 0 < t(k) < 1 is 
the plasmonic spectral weight, and f l(k) is the plasmonic dis¬ 


persion with fl(0) = ar s / (37t) ep, where cf is the Fermi 
energy. For the numerical integration we define the bare time- 
ordered Green’s function as 26 


G 0 -(fe, w) 


B(k) | A(k) 
co — ek — irj to — €k + iij 


(58) 


In non-interacting systems B{k) = rik and A{k) = 1 — Uk 
with rife denoting the occupation of the state with momentum 
k. The frequency integrations can be done completely ana¬ 
lytically (facilitated by the MATHEMATICA computer algebra 
system (CAS)) whereas for the remaining momentum inte¬ 
grations one has to rely on numerics 38 . The starting point are 
the four diagrams depicted at Fig. 15. The momentum flows 
are explicitly shown. There are in general many possibilities 
to assign momenta to propagators. Our choice is dictated by 
the matter of convenience and is not unique. In order to fur¬ 
ther simplify notations we adopt the following short forms: 
A, = A(xi), Bi = B(xi), Ci = \ v(yi)w(yi ), and introduce 
the function 


n t {a,n) 



B j 

CL T fl — Ci 


We also recall that |w|) = B < (k, — |w|). Therefore, it 

is sufficient to consider only one case, e.g., co > 0. 












14 


The results of frequency integration are 

' d 3 x\ 


^a{z,0=-nj 

&b( z ’0=-nj 
B?( z ,0=-nj 
(^C)=-7Tj 


(2tt) 3 


A^B^S^C, — £4 + £ 3 ), 


'd 3 X\ I" d 3 yi Tieie 4, ^1) - ?4(e3, ^1) 
(2n) 3 J(2nY 


e 3 — e 4 — e 5 + e 6 


' d 3 x\ 

WP. 

' d 3 x 1 


'd 3 t /2 di 2 {t^ ^ 2 ) — Hi(£ 3 , ^ 2 ) 


(2tt) 3 


— £4 + £ 3 ), 


C , 2^3-64 < 5(C — e 4 + £ 3 )) 


£3 - £4 - £l + £2 
’ d 3 2/2 Beici, f2i) — %5(e3, ^ 1 ) T~L 2 (e4, ^ 2 ) — %i(e3) ^ 2 ) 


(2n) 3 J (2ir) 3 J (27r) 3 £3 - £4 - £5 + 


£6 


£3 - £4 - £1 + £2 


(59a) 

(59b) 

(59c) 

CiC'2A 3 B4 l 5(C-£4 + e3)(59d) 


They are quite general and can be used to obtain e.g. plasmonic contribution. If, however, results for bare Coulomb interacting 
lines are needed we take the limits and obtain: 


Ba( Z ,0=-Kj 

^(^C)=-7tj 

®d 0>C)=-7tj 


" d 3 xi 
(2^)3 

' d 3 Xi 

(2^)3 

" d 3 X\ 

(2tt)3 
" d 3 Xi 


^-3^4^(C — e 4 + £ 3)1 

f d 3 yi , N ^6 A 5 
(2^)3 

f d 3 y 2 


(2tt) 3 

^ 3 yi 


^(yi) 

^(2/2) 


£3 — £4 — £5 + e 6 

A2 — -Ai 


£3 - £4 - £l + £2 


A 3 .B 4 < 5 ((( — £4 + £ 3 ), 
A^,B^3{C, — £4 + £ 3 ), 


'rf 3 y 2 


(27t) 3 J (27t) 3 J (2ny 


v{y 1 )v{y 2 ) 


Ak — Ac, 


A 2 — Ai 


£3 - £4 - £5 + £6 £3 - £4 - £l + £2 


(60a) 

(60b) 

(60c) 

-A 3 £?4<5(C-£4 + £3)- (60d) 


Notice that the spectral functions in Eqs. (59) and (60) are 
denoted by the same symbol because their type can always 
be inferred from the context. In these equations A, B, and 
£ quantities are labeled by the momenta as shown at Fig. 15. 
For instance 

f 1 \x 1 -y 1 -y 2 \< k F , 

5 \ 0 jxi - y 1 - y 2 j > k F , 

and £5 = (a?i — 2/1 — y 2 ) 2 / 2. B a (z,0 = B (0) {z, () is 

obviously the Findhard polarization function. Bb{z,C,) and 
B c (z, 0 differ only by the permutation of indices and, there¬ 
fore, are equal in view of the left-right symmetry of the cor¬ 
responding diagrams. There are no other topologically iden¬ 
tical diagrams of the first order, i.e. Bb(z,Q + B c (z,Q = 
B^\z,Q. There are two more partitions of the second order 
having the same topology as Bd(z,Q■ They have different 
combinations of pluses and minuses assigned to the vertices 
and are not considered here, hence Bd{z, £) = £). 

Before analyzing numerical results let us first notice a dif¬ 
ferent proportionality of each perturbative order to the elec¬ 
tron density. Because we scaled all spectral functions such 
that £r 0) is density independent it is easy to see that B 1 '= 
0((ar s ) n ). Thus, the sum of four terms in Eqs. (60) is given 
by a quadratic polynomial in terms of ar s . The requirement 
of its positivity leads, therefore, to the following inequality 

V= (B<(z,0) 2 -B<(z,0B<(z,0 < 0 , ( 61 ) 

which should hold for all values of frequency and momenta 
for at least one density. This ensures the positivity of the spec¬ 
tral function at all densities. Mathematically, inequality (61) 


is the Cauchy-Schwarz inequality applied to the integrals of 
half-diagrams. 

The 1st and 2nd order expressions (Eqs. (60)) as well as 
the analytic result of Holas et al. [17] (Appendix A) suf¬ 
fer from the logarithmic singularities of the integrated func¬ 
tions. This poses additional challenges for numerics. We 
tackle this problem by introducing a small imaginary part 
iS into the energy denominators. The Monte-Carlo integra¬ 
tion is performed as in Ref. 26 with the use of Mersenne 
twister 19937 random number generator 39 . Additional com¬ 
plication arises due to the use of bare rather than the plasmon- 
ically screened Coulomb interaction, i.e. even large momen¬ 
tum transfers are possible. In this case we can still map the 
integration variables to the [0,1] interval by the logarithmic 
scaling. For instance, we represent the vector yi as follows 
ri = — plog(xi), cos(#i) = —1 + 2^2, ft 1 = 277x3, where 
p is some suitably chosen constant 40 and Xi € [0,1]. Thus, 
d 3 yi = 4^dxidx 2 dx 3 . 

Fet us look now at the results of the Monte-Carlo integra¬ 
tion for a density of r s = 3. The full momentum (0.1 k F < 
k < 2.5 k F ) and energy (0 < u> < 12e F ) resolved spectral 
functions are shown in Fig. 16. On this graph the numeri¬ 
cally produced points are projected on the analytically known 
results depicted as surfaces. For the second order analytical 
expressions are not known, therefore, only points are shown. 
At small momenta the spectral functions diverge, therefore 
we introduced some truncation. The sum of all tree terms, 
Bpsd, is always positive as an example of a cross-section at 
k = 1.2 k F in Fig. 17 demonstrates. For some frequencies, 
however, Bpsr>{k, x 1 ) is rather small and lies within the error 
bar of the Monte Carlo simulation. 

The domain where B^\k, ui) is negative is not bounded. 
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FIG. 16. (Color online) The scaled spectral function for the polariz¬ 
ability V of homogeneous electron gas at r 3 = 3. Surfaces denote 
exact analytical results: Oth order is given by the Lindhard function 
(Eq. (56)) and 1st order is computed according to Id integral repre¬ 
sentation of Holas et a/. 17 . Dots denote numerical results obtained 
by the Monte-Carlo calculation. 



FIG. 17. (Color online) Cross-sections of the results in Fig. 16 for 
the momentum value of k = 1.2kF and density r a = 3. Top left: 
zeroth order contribution £^ ()) . Top right: first order contribution 
S (1) . Bottom left: second order contribution S (2) . Bottom right: the 
sum of all three contributions Bpso which is positive. Dots (blue) 
denote numerical Monte Carlo results. Solid lines (red) stand for 
analytical results. 


see again Fig. 14. Thus, corrections originating from 
{k,oj) qualitatively modify the behavior of the spectral 
function at any momentum. This, in view of the Hilbert trans¬ 
form in Eq. (17), leads to a modification of the real part of the 
response function. The general conclusion is that the cutting 
procedure for PSD spectra works as it should. The addition of 
the second order vertex diagram correctly removes the nega¬ 
tive parts of the response function to first order in the vertex. 


Cancellation between vertex corrections and self-energy 
insertions 

Although our example is the simplest one that illustrates 
the PSD diagrammatic theory, it is too simple from a physical 
point of view. The main reason is that we used bare propa¬ 
gators Go and bare interactions v. The first order vertex di¬ 
agram is not the only first order diagram as we missed two 
first order diagrams that contain exchange self-energy inser¬ 
tions. If we had used an expansion in dressed Green’s func¬ 
tions G and dressed interactions W such diagrams would not 
appear as in that case we could restrict ourselves to skeleton 
diagrams. However, in an expansion in Go and v they be¬ 
come relevant. In particular the self-energy diagrams lead to 
a cancellation of the divergent small fc-behavior of the spec¬ 
tral function For the case ui = 0 this as been ex¬ 

plicitly demonstrated by Engel and Vosko 2 ’. It would there¬ 
fore be a natural step to include the first order self-energy 
diagrams as it would, for instance, guarantee the existence 
of a gradient expansion for the exchange-correlation energy, 
see below. Taken together these diagrams can be expanded 
in a power series in terms of the momenta which plays an 
important role in determining the gradient expansion of the 
exchange-correlation energy functional in density functional 
theory. It is well known that in the lowest order the gradi- 
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ent correction to the exchange-correlation functional can be 
expressed as 23,41 

AE„[n 0 ,fa(q)] = \j 

where 5n is the density variation with respect to the density 
no of the homogeneous system, and A"P(q, 0) denotes cor¬ 
rections to the response function from the first and higher or¬ 
der diagrams. Therefore, the inclusion of B r2> will certainly 
modify the gradient expansion coefficients, e.g, due to En¬ 
gel and Vosko 2 '. However, to have a well-defined gradient 
expansion one has to add the self-energy diagrams as well. 
Unfortunately, a simple addition would destroy the positiv¬ 
ity of the resulting spectral function again. This was noticed 
e.g. by Brosens and Devreese 22 and is illustrated in Fig. 18, 
where we also display the contribution of the self-energy dia¬ 
gram (B^ 1,Se \k, w)) calculated using the analytic expression 
of Holas et al. 11 . Therefore, if we desire to include the self¬ 
energy diagrams and still wish to guarantee positivity we have 
to apply our PSD theory and consider an extended set of di¬ 
agrams. The minimal set that achieves this goal is displayed 
in Fig. 19 which apart from additional self-energy diagrams 
also contains mixed self-energy and vertex diagrams. Rather 
than developing codes to evaluate these additional diagrams 
we found it more worthwhile to explore approximations that 
involve dressed Green’s functions and interactions. First of 
all, the dressing of the interaction reduces the singular behav¬ 
ior of the diagrams and secondly they reduce the number of 
diagrams to be evaluated since we can then stick to skeletonic 
diagrams. However, since this requires an extensive discus¬ 
sion by itself we will address this topic in a future publica¬ 
tion. An alternative route was undertaken by Brosens, De¬ 
vreese and Femmens in a series of works 18 - 19 - 42,43 using the 



FIG. 18. (Color online) Spectral functions at momentum k = 2.5 kp 
and density r s = 3. (a) In addition to previously considered ze¬ 
roth order (S (0 \ dotted) and first order (Z3 (1 \ short dash) a contribu¬ 
tion of the diagrams with self-energy insertion (£! (1 ’ Se ' ) , long dash) is 
shown, (b) Sum of three terms in panel (a), (c) Second order contri¬ 
bution, B 1 ' 2 ' 1 . (d) Bpsd including the first order self-energy diagrams 
(dots), solid line as in (b), shaded area denotes second order contribu¬ 
tion. In order to cancel small negative spectral function at the edge of 
particle-hole continuum (u> = 1.25 £f) inclusion of more diagrams 
as shown at Fig. 19 is required. 



FIG. 19. (Color online) "Ppsd obtained from the cutting procedure 
including the first order self-energy diagrams. 


variational solution of the linearized equation of motion for 
the electron density distribution function. They have shown 
that the first term in a series expansion of the variational re¬ 
sult for the local field factor yields the lowest order diagram¬ 
matic result. Since such solution also contains higher order 
terms the spectral function is positive. However, it is not free 
from singularities suggesting again possible benefits of work¬ 
ing with dressed Green’s functions and interactions. 


VII. CONCLUSIONS AND OUTLOOK 

Vertex corrections in diagrammatic approximations to the 
polarizability are known to be crucial for capturing double and 
higher particle-hole excitations, excitons, multiple plasmon 
excitations, etc. as well as for estimating excitation life-times. 
However, the straighforward inclusion of MBPT vertex dia¬ 
grams can lead to negative spectra, a drawback which discour¬ 
aged the scientific community to develop numerical recipes 
and tools for the evaluation of these diagrams in molecules 
and solids. In this work we provided a simple set of rules 
to select special combinations of diagrams yielding a positive 
spectrum. 

In our formulation every MBPT diagram is written as the 
sum of partitions, and every partition is cut into half-diagrams. 
We recognized that they are the half-diagrams the fundamen¬ 
tal quantities for a PSD expansion. In fact, the sum of squares 
of half-diagrams corresponds to a special selection of parti¬ 
tions which is PSD by construction. The requirement of pos¬ 
itivity on the spectrum is important not only for the physi¬ 
cal interpretation of the results but also for the correct an¬ 
alytic structure of the polarizability. We demonstrated that 
a PSD polarizability cannot generate a (retarded) density re¬ 
sponse function with poles in the upper-half of the complex 
frequency plane. This is a critical property to converge self- 
consistent numerical schemes. Although the PSD diagram¬ 
matic expansion put foward in this work applies equally well 
to bare as well as dressed Green’s functions a word of caution 
is due in the dressed case. The gluing of skeletonic, i.e., self¬ 
energy insertion-free, half-diagrams can lead to nonskeletonic 
polarizability diagrams. In order to avoid the double count¬ 
ing of some of the diagrams it is therefore necessary to use 
Green’s functions dressed with self-energy diagrams distinct 
from those appearing in V. 


















A natural way to sum to infinite-order a subclass of polariz¬ 
ability diagrams is through the BSE, an integral equation with 
kernel given by the functional derivative of the self-energy 
with respect to the Green’s function. Due to the popularity of 
the BSE we also addressed the issue whether the polarizability 
which solves the BSE is PSD for conserving self-energies, and 
found a negative answer. The counter example is provided by 
the self-energy in the second-Born or GW approximations. 
Noteworthy these self-energies yield a positive spectrum for 
the Green’s function 26 ; therefore neither a conserving nor a 
PSD self-energy does necessarily generate a PSD polarizabil¬ 
ity through the BSE. 

The simplest approximation with vertex corrections is the 
first-order ladder diagram. This diagram has been calculated 
both in finite and bulk systems and it is known to be not 
PSD. How to include vertex corrections without altering the 
positivity of the spectrum has been a long-standing problem, 
which we have solved in this work. By adding a partition 
of the second-order ladder diagram we obtained the simplest 
PSD approximation with vertex corrections. We then evalu¬ 
ated this approximation in the 3D homogeneous electron gas 
and confirmed numerically the correctness of our PSD theory. 
We stress again that the PSD property alone does not neces¬ 
sarily guarantee physically meaningful spectra. In fact, the 
PSD spectrum with vertex corrections has un unphysical di¬ 
vergency at zero frequency and momentum. The inclusion of 
bubble diagrams with first-order exchange self-energy inser¬ 
tions removes this divergency but destroys the PSD property. 
We worked out the minimal set of diagrams to turn this ex¬ 
tended approximation into a PSD one. The evaluation of the 
resulting extra diagrams is within reach of our code but it re¬ 
quires a considerable numerical effort and it goes beyond the 
scope of the present work. 


VIII. ACKNOWLEDGMENTS 

AMU would like to thank the Alfred Kordelin Foundation 
for support. GS acknowledges funding by MIUR FIRB Grant 
No. RBFR12SW0J. YP acknowledges support by the DFG 
through SFB762. RvL would like to thank the Academy of 
Finland for support. 


Appendix A: Numerical evaluations of the analytical 
expressions for the first order polarizability 


As in the rest of the text we measure the momentum and 
the energy in units of the Fermi momentum fcp = l/(ar s ) 
and the Fermi energy cf = kp/2. In older papers k‘ 2 F as 
the energy unit was used 17 . This must be taken into account 
when comparing. Also notice that Engel and Vosko measured 
momentum in terms of 2k f- It is natural because the first 
order polarizability has a logarithmic singularity at this point. 
The imaginary part of the dielectric function resulting from 
the first order polarizability is given by 


Ime((7, w) 


8 a 2 r 2 

7T fc 4 


(F E -(t,|) + r-(t,|)), (ad 


17 


whereas the scaled spectral functions considered in Sec. VI 
read 

( A2) 

= <a3) 

with the function F non-zero at two domains in k - oj plane 
(restricted by the Heaviside ^-functions): 


F Ex > Se (fc,w) =0 





k 

2 


(j)Ex,Se 


*"*(!-i*) 
(-H4 <«> 


and in turn 

<b s >, k ) = kf L ([(k + v) 2 + (1 - v 2 ) ] 1/2 ) 

~{k + v)f L {k + v) + vf L (v), (A5) 
$ Ex (+ k) = -Gi(zz) + G 2 (v + k, 1 - i/ 2 ), (A6) 

where /l is the Lindhard function 


1 1 — z 2 

■«--> = 2 + T^‘» g 


2 + 1 


1 - 2 


and G\ is (cf. Eq. (2.15) of Holas et a/.in Ref. [17]) 

G.M = j(l -+s(l4) - + (A7) 

x ((1 - v) log(l - v) + (1 + v) log(l + v) — 2 log 2), 

with g(z) = Li2 (—2) — Li2( — I/2) represented in terms 
of the polylogarithm functions. The second function is more 
involved, it is given by the Hilbert transform which has to be 
computed numerically: 


G 2 {x,y) 


1 

4 



T(£,x,y) 

Z-x 


(A8) 


If | a; | < 1 the simplest way to avoid singularity is to exclude 
a small (\x a — Xb\ < e) interval x £ (, x a , Xb) C (—1,1) from 
the integration. Finally, the T{^,x,y) function is defined as 
(cf. Eqs. (2.18-2.23) of Holas et al. in Ref. [17]): 

T (£, x ,y) = (t| - y + (! - £ 2 ) (1 - *)) 

x log (2t (l — £ 2 ) + ai) — y log (ai) 

+ (1 - e) (t - log(t) - 1) (A9) 


in terms of auxiliary functions (7 = x 2 

a 2 = 4(1 - £ 2 )(£ - x) 2 , a 3 = 2£(£ - 

and 


+ y,ai= 2(£ - a:) 2 , 
- - 1 ’ A = (7+a 3 ) a) 


t(€,x,y) 


(t-x ) 2 

03+7 

03+7 yi+A —1 

(l-? 2 ) 2 


|A| < e, 
|A| > e. 


(A10) 


From the imaginary part of the polarization function the real 
part can be computed through the Hilbert transform. For the 
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static case co = 0 we have: 


where for a(q) and b(q) there are analytic expressions due to 
Engel and Vosko 2 ’: 


Rex (1) (fc,0)= 


2 

7r 3 fc 2 
1 
7T 3 


+ (fc+2) dio 

CO 
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a I - 


F E *(k,“)+F s °(k,“) 
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(All) 
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log 2 


1 + X 


■log 


1 — X 

+ 9 


dx 


1-9 


-6(9). 


(A12) 


dx 


(A13) 


Numerically Eq. (A13) is much faster than the Hilbert trans¬ 
form (A11). However, it is good to know that both ways yield 
identical results that also agree with our Monte-Carlo simula¬ 
tions. 
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